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ABSTRACT 

Using a version of the ZEUS code, we carry out two-dimensional simulations 
of self-gravitating shearing sheets, with application to QSO accretion disks at a 
few thousand Schwarzschild radii, corresponding to a few hundredths of a parsec 
for a 100-million-solar-mass black hole. Radiation pressure and optically thick 
radiative cooling are implemented via vertical averages. We determine dimen- 
sionless versions of the maximum surface density, accretion rate, and effective 
viscosity that can be sustained by density-wave turbulence without fragmenta- 
tion. Where fragments do form, we study the final masses that result. The 
maximum Shakura-Sunyaev viscosity parameter is approximately 0.4. Fragmen- 
tation occurs when the cooling time is less than about twice the shearing time, 
as found by Gammie and others, but can also occur at very long cooling times in 
sheets that are strongly radiation-pressure dominated. For accretion at the Ed- 
dington rate onto a 10 8 solar-mass black hole, fragmentation occurs beyond four 
thousand Schwarzschild radii, r§. Near this radius, initial fragment masses are 
several hundred suns, consistent with estimates from linear stability; final masses 
after merging increase with the size of the sheet, reaching several thousand suns 
in our largest simulations. With increasing black-hole mass at a fixed Eddington 
ratio, self-gravity prevails to smaller multiples of rg, where radiation pressure 
is more important and the cooling time is longer compared to the dynamical 
time; nevertheless, fragmentation can occur and produces larger initial fragment 
masses. Because the internal thermal and gravitational energies of these massive, 
radiation-pressure-dominated fragments nearly cancel, small errors in energy con- 
servation can cause spurious results such as spontaneous dissolution of isolated 
bodies, unless special care is taken. This is likely to be a challenge for all eulerian 
codes in self-gravitating regimes where radiation pressure dominates. 

Subject headings: accretion, accretion disks — galaxies: active — gravitation — 
methods: numerical — stars: formation 
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Introduction 



Active galactic nuclei (AGN) are powered by accretion onto a supermassive black hole 
(SMBH). Probably but less certainly, accretion occurs via gaseous disks, and furthermore via 
a local balance between inward advection of angular momentum and outward transport by 
torques due to magnetohydrodynamic turbulence, spiral density waves, or magnetized winds. 
Though mostly indirect, the evidence for disks at distances < 1 pc from the SMBH is ex- 
tensive. Theoretically, disks naturally result from the combination of energy dissipation and 
angular-momentum conservation; also, disk accretion down to the marginally stable orbit 
naturally accounts for conversion of ~ 10% of accreted mass to radiation, as req uired by com 



parisons of the integra ted light from AGN with the integrated mass in SMBHs (jSoltanlll982 



Yu fc Tremaind |2002| ). Phenomenologically, disk accretion is consisten t with the evidence 



for an isotropic emission as embodied in the unification model for AGN f lMiller fc Antonucci 



19831 ). In a few highly selected AGN, VLBI observations of maser emission strongly indicate 
the presence of a gaseous disk in keple rian rotation at a few tenths of parsecs from the SMBH 
jMivoshi et al.lll995l : IKuo et al.lfioioh . 



A longstanding difficulty with the disk paradigm for AGN is to explain how gas is trans 



porte d from galactic to sub-parsec scales without turning entirely into stars (IShlosman et al. 



19901 ). Disk accretion in local thermal equilibrium is prone to self-gravity because the ver- 
tical thickness tends to be small and the inflow speed slow, so that the density of the gas 
must be large to support inferred mass accretion rates. On scales > GMbh / cr circ,gai ~ 10 P c ' 
external gravitational torques due to merging galactic nuclei or stellar bars may speed the 
inflow. On scales < 10 3 GM BH /c 2 ~ 10~ 2 pc, the tidal field of the SMBH is large enough and 
the temperature of a standard viscous disk is high enough so that self-gravity is typically 
slight. This leaves a broad range of radii, typically 0.01-10 pc or 10 3 -10 6 Schwarzschild radii 
rs, over which a standard viscous disk in steady-state accretion would b e self-gravitating , 
with a Toomre Q parameter that decreases rapidly with increasing radius (|Goodmanll2003l ) . 
Gravitational fragmentation, leading to star formation within the disk, appears to be a nat- 
ural outcome. Inde ed, stellar disks a re common within a f ew parsecs of the S MBH in nearby 
quiescent galaxies ( ILauer et al.ll2005l ) and active Seyferts f IDavies et al.ll2007l ). w hile our own 
Gala c tic Center contains what appear to be main-sequence B stars at ~ 10 -2 pc (IGhez et al. 
20031 : iMartins et alJliooi l. 



There is a general belief t hat star formation and accretion somehow regulate and sup- 
port one another in AGN (e.g.. ICollin fc Zahnlll999l ; iThompson et al.ll2005l ). But, we believe, 
a convincing model of how this works is still lacking. Beyond the usual perplexities attend- 
ing "normal" star-formation, sub-parsec AGN disks pose severe problems of energetics and 
stability. Specific orbital energies are not small compared to those released by stellar evolu- 
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tion (~ 10~ 3 c 2 ), while orbital and thermal timescales are much shorter than the minimum 
main-sequence lifetime (~ 10 6 yr); both comparisons would appear to make a stable feed- 
back between gravitational collapse and stellar energy inputs more difficult than in giant 
molecular clouds. Also, it is not clear that AGN disks should form stars of normal mass. 
Evidence exists f or a somewhat top-heavy stel l ar mass function am o ng; the young stars in the 
Galactic Center ( iNayakshin &: Sunyaevl l2005t iBartko et al.ll2010l ). iGoodman fc Tan! (120041 ) 
argued that objects of up to ~ 10 5 M might form, this being the so-called "isolation mass" 
(a dynamical scale borrowed from theories of planet formation) appropriate to radii ~ 10 3 r$ 
in a typical bright AGN. 

In view of the physical and observational difficulties of this subject, a complete theory 
may not be available for many years, but progress has been made in unders tanding the role 
of self-gravity in promoting accretion, and the conditions for fragmentation. iGammiei (120011 ) 
showed via idealized two-dimensional simulations that disks (or rather, shearing sheets) 
subje ct to coo l ing ca n maintain themselves in a state of marginal linear stability according 
to the iToomrd (119641 ) criterion provided that the product Qt c cooling time and orbital angular 
velocity is somewhat greater than unity, where t c is the timescale on which the gas would 
radiate its thermal energy in the absence of heating processes. Gammie's models with 
Qt c > 1 reached a statistical steady state with nonlinear density waves and shocks sufficient 
to offset the cooling. Since the energy source that supports ongoing mechanical dissipation 
is ultimately differential rotation, the wave turbulence must transport angular momentum 
outward, and since Gammie's models were effectively local this transport can be described 
by a Shakura-Sunyaev viscosity parameter a ~ (f2t c ) _1 , the exact value depending upon the 
assumed ratio 7 of specific heats of the gas. Below a critical value flt c ~ 3, Gammie's disk 
fragmented into gravitationally bound objects. His conclusions have been confirmed, with 
some variations i n the critical values of Q t r and a, by subsequent simulations with more 
realistic cooling ([Johnson &: Gammid 120031) and with global, three-dimensional geometries 
faice et alJl200a lLodato fc Ricell2004l ). 

Recently, ( Hopkins fc Quataert 2010bl lal hereafter HQ) have argued, based on SPH 
simulations buttressed by analytic arguments, that global, nonlinear, non-axisymmetric den- 
sity waves accompanied by shocks can support accretion rates as high as ~ 1OM yr -1 at 
< 0.1 pc. One-armed spirals (azimuthal mode number m — 1) are particularly prominent 
in their simul ations, probably becau se such modes cohere most easily in near-keplerian po- 
tentials (e.g., iLee fc Goodmanlll999l and references therein). It is true that nonlinear global 
spirals can in principle exert much larger torques on the gas than is possible with a lo- 
cal effective viscosity, by a factor ~ {rjj%}a~\ where h is the disk thickness and a is the 



Shakura-Sunyaev viscosity parameter (lGoodmanll2003l . hereafter G03). HQ's results sidestep 
rather than solve the problem of local self-gravity, however. Their simulations employ a su- 
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perthermal effective sound speed that is intended to represent unresolved turbulence, and 
which is large enough to suppress local instabil ity. This device was originally developed to 
parametrize stellar feedback on galactic scales (jSpringel et al.ll2005l ). but for the reasons of 
energetics and timescales mentioned above, one may question its applicability to the scales 
of interest to us, 0.01-0.1 pc from the SMBH. In any case, global waves, particularly HQ's 
near-stationary m = 1 waves, could not be observed in shearing-sheet simulations such as 
those of our paper. 



In all the simulations following iGammid (120011 ) 's original work, one issue that has not 
been directly addressed is the role of radiation pressure in the equation of state, perhaps 
because most of the applications have been to protostellar disks or to the Galactic Center. 
For accretion rates and black-hole masses characteristic of bright AGN, however, radiation 
pressure is still important at the minimum radius where self-gravity sets in (G03). As is well 
known, the ratio of radiation pressure to gas pressure is a monotonically increasing function 
of mass for optically thick, nondegenerate bound objects in hydrostatic equilibrium; and 
the effective adiabatic index Ti = {d Inp/d In p)$ a correspondingly decreasing function. 
This suggests that as objects gain mass through accretion or merging in a fra gmenting 
AGN disk, they become increasingly suscep tible to rap i d coll apse. As noted by iGammie 
(120011 ) and confirmed in three dimensions by iRice et al.l (120051 ) , the critical value of Qt c for 
fragmentation depends upon the effective equation of state relating pressure (p) to mass 
density (p), or height-integrated pressure, P, to surface mass density, S. In particular, if 
I2D — (dhxP/dki S)g < 3/2, then collapse may occur without any cooling (Qt c = 00). This 
perhaps is why Gammie chose to do his simulations with / -f 2 D — 2. The correspondence 
between T± and 72D depends upon how the vertical thickness of the disk responds to changes 
in surface density: if the response is hydrostatic and strongly self-gravitating, then 72D = 3/2 
corresponds to Ti = 4/3, the value for a spherical body supported entirely by radiation 
pressure. 

Most of these points reg arding the importance of radiation pressure were made by G03 
and lGoodman Sz Tan! (120041 ) via analytical arguments; the principal goal of the present work 
is to illustrate them through numerical simulations. We originally hoped also to demon- 
strate growth up to the isolation mass (appropriately redefined for the shearing sheet), but 
while we do demonstrate growth well beyond the mass scale associated with linear insta- 
bility ("Toomre mass"), numerical difficulties described below prevented us from making a 
systematic study of the ultimate masses as a function of shearing-sheet control parameters. 



The structure of this paper is as follows. In §2.11 we present the adopted equation of 
state for our height-integrated, two-dimensional calculations; some mathematical details are 
deferred to an Appendix. Our equation of state incorporates both radiation pressure and self- 
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gravity under the assumption of local vertical hydrostatic equilibrium. In §2.2[ we review the 
basic equations for the shearing sheet. Our cooling prescription, which is based on a simple 
algebraic approximation to vertical radiative transfer, is presented in §2.31 We describe our 
computational units and control parameters in §2.4} and our numerical methods in §EJ In 
§HJ we present results from representative simulations on both sides of the fragmentation 
boundary. A general picture of the disk based on our simulations is described in §4.31 A 
summary and discussion of our main results, and some speculations concerning observable 
consequences, are given in §SJ 



This section presents the height-integrated physical model upon which our numerical 
simulations are based. 



Let (r, <p, z) be cylindrical coordinates such that the central black hole, with Schwarzschild 
radius r s = 2GMbu/c 2 , lies at r = 0, and the disk midplane at z = 0. At radii r ~ 10 3 - 
10 4 rs, the disk is expected to be quite thin, with a half thickness h ~ lCT 2 r (G03). The z 
coordinate is therefore eliminated from the governing equations of our numerical simulations 
by vertical integration so that, for example, pressure p and mass density p are replaced by 
height-integrated pressure P and surface density E: 



By an "equation of state," we mean a functional relation among P, S, and other thermo- 
dynamic variables. To obtain such a relation, we make a number of simplifying assumptions 
about the thermodynamics of the gas and about the distribution of p and p with z. Since the 
regions of the disk with which we are concerned are probably very optically thick, the gas 
and radiation temperature are taken equal. Vertical hydrostatic equilibrium is assumed, and 
magnetic and turbulent contributions to the zz component of the stress tensor are neglected, 
so that the disk thickness is supported entirely by the sum of gas and radiation pressure, 
p = p r + P g - The gas pressure fraction 



2. A shearing-sheet model for AGN disks 



2.1. The equation of state 





(1) 



Pg+Pr 



(2) 
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is presumed constant with height but allowed to vary with (r,(f),t). In combination with 
vertical hydrostatic equilibrium, constant /3 implies constant nF z /g z , where k is the opacity; 
g z is the vertical component of the gravitational field; and F z is the vertical component of the 
radiative heat flux. The vertical runs of density and pressure are then related by a polytropic 
relation, 



P 



3 / k B \ 4 1 -0 



1/3 



(3) 



a \jim p ) f3 A 

with molecular weight p ~ 0.62, as appropriate for a fully ionized gas of near-solar metallicity. 

Once P is specified, P can be found in terms of S by inserting eq. ([3]) into the equation of 
vertical hydrostatic equilibrium. But when self-gravity is included [via the one- dimensional 
form flAip of Poisson's equation], the relationship that results is implicit. We relegate the 
mathematical details to the Appendix and simply quote the main results here. 

It is convenient to introduce the quantity 

Q 2 

where p(0), which is shorthand for p(0,r,(j),t), is the mass density at the midplane, and 

f2 = (GMbh/^ 3 ) 1 ^ 2 is the orbital angular velocity. We use this notation because Q as we 

have defined it is usually numerically comparable to Toomre's stability parameter Qt = 

VlCs/itGYj, to which it would reduce if the effective thickness S/p(0) of the disk were given 

i li 

by 2C S /Q in terms of the effective horizontal sound speed C s = (dP/dT,)^ . The actual 
value of the effective thickness is somewhat different from 2C S /Q, however, not only because 

1 /2 

of vertical variations in the true three-dimensional sound speed c s = (dp/dp) s , but also 
because of the self-gravity of the disk. The quantity (jl]) is more directly related to the Roche 
criterion for an object of mean density ~ p(0). 

It is shown in the Appendix that 

_ ^Qh{Q) G 2 S 3 

16[/ 3 (Q)] 3 fi 2 ' U 

(ttG) 7 / 3 Q 4 / 3 £ 2 



and 



K V) 2 "/3 08/3 lh{Q) f • (6 ) 

(7 ) 



3ir 2 h(Q)QG 2 Z 
in which 

/ (O) ~ °' 323 / (O) ~ °' 287 
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Equation (JSJ) determines P in terms of E and Q, but Q is not a conserved or primitive 
variable in our dynamical equations. Instead, we have the thermodynamic internal energy 
per unit area, U. Since the internal energy per unit volume u = (3/2)p g + 3p r = 3(1 — (3/2)p 
under our assumptions of complete ionization and equal gas and radiation temperatures, and 
since j3 is assumed independent of z, it follows that 

ff-(i-|)p. (9) 

Since K{0) is the function ([3]), equations (jHJ) and (J7|) determine (3 and Q in terms of E, [/, 
and fi, and then eq. OH]) gives P. It can be shown from eqs. ©, 0, and © that P oc E 3 / 2 
when both /J and Q are ^ 1; i- e - the effective 2D adiabatic index approaches the critical 
value of 3/2 at which nonrotating bound fragments can collapse indefinitely without cooling. 

For the calculation of the local cooling time, we sometimes require the physical density 
and temperature since we include a Kramers component in our opacity law ( §2.2p . The 
mid-plane density p(0) follows from Q via eq. 01]), and the mid-plane temperature T(0) is 



T(0) 



3k B 1-/3 ' 

a 

ajj,m p p 



1/3 



(10) 



2.2. Dynamical equations 



Following iGammid (120011 ) . we describe the local dynamics of the disk in a shearing sheet 
approximation, in which x = r — r and y = r [<f) — fi(r )t] are pseudo-Cartesian coordinates 
centered on a circular orbit of radius r$. The equations of motion are 



0. 



Dv 
~Dt 



VP 



2fie z x v + SQ 2 ^ - V$, 
A. 



-PV ■ v 

V 2 $ = 47rGE(5(^). 



(11) 

(12) 

(13) 
(14) 



The cooling function A represents radiative losses from the surface of the disk, as described 
in §2.31 If A = 0, these equations have steady solutions in which E and $ are constants and 
v = — ^Qxe y . 



A useful diagnostic quantity is the potential vorticity 

_ V x v + 2f2 



(15) 
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This is conserved following the fluid, D^/Dt = 0, in the absence of shocks, viscosity, or 
uneven cooling, so that the equation of state is effectively barotropic, P = -P(E). 

We do not include any explicit viscous terms. Since we do not resolve the vertical 
dimension of the disk/sheet, we could not represent magnetororational instabilities (MRI) 
directly but would have to parametrize the magnetic stresses in terms of £ and U ; such 
parame trizations are prone to thermal and viscous inst abilities where radiation pressure dom- 
inates ( iLightman fc Eardleyill974t iHirose et al.ll2009l ). Furthermore, the effective Shakura- 
Sunyaev parameter provided by density waves and shocks in our simulations is usually larger 
than omri ~ 10~ 2 . There is, howe ver, an artificial viscosity included in ZEUS to mediate 
shocks (e.g.. IStone fc Normanlll992[ ). 



2.3. Cooling Fuction 



Following I Johnson fc Gammid (120031 ). the cooling function A in equation (TT3"j) represents 
radiation losses from the surface of the disk. Since we do not resolve the disk thickness we 
adopt a standard algebraic approximation for the vertical radiative transfer: 



A 



laTf 



16 



cff 



aT 4 (0) r + 



1 



T 



-1 



in which 



1 



T 



£K[p(0),T(0)]. 



(16) 



(17) 



approximates the local optical depth from the midplane to the surface. Usually r > 1 on 
average in those parts of AGN disks that we wish to model, but in case fragmentation should 
lead to gaps in the sheet, the form of equation 0161) is ch osen so that A oc r in the optically 



thin limit (e.g. Hubeny 1990 ; Johnson fc Gammiel 



2003|). 



At radii 10 3 -10 4 rs, the mid-plane temperature T(0) is typically 10 4 -10 5 K for near- 
Eddington accretion rates (GO03). If we assume Q = 1, the mid-plane density will be 
p(0) ~ 10 _9 gcm -3 . In this density and temperature range, the dominant opacity is electron 
scattering, which is almost constant, justifying the factor 1/2 in equation (I17p . We include 
a Kramers opacity, however, which often begins to be important beyond ~ 5000rg in our 
simulations. An analytic approximation to the opacity sufficient for our purposes is therefore 



+ k k = 0.2(1 + X) + 4 x 10^° (1 + X)(Z + 0.001) 



P 

^3.5 



(18) 



all quantities being evaluated in in cgs units (B. Paczynski, private communication). The 
mass fractions of hydrogen and metals are taken at their solar values, X = 0.7, Z = 0.02. 
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There is evidence that the broad-line gas in quasar stell ar objects is more metal-rich than the 



sun (e.g. lHamann Sz Ferlandlll993l ; iDhanda et al.ll2007l ); however, this would not much affect 
our results because electron scattering opacity dominates as long as the metallicity of the 
broad-line gas is smaller than ~ 10 of the solar value, at least for near-Eddington accretion 
onto black holes of masses > 10 8 M Q in the range of radii where self-gravity is important 
but fragmentation is avoided. Kramer's opacity gains in importance with increasing radius, 
decreasing M B h, and decreasing M. 



The local cooling time is defined as 



U 



(19) 



For Q > 1, the half-thickness h ~ y'P/E. Tt follows that t c ~ kT 4 (0)/(c^) 2 for < 1, and 
that t c oc kS 2 /T 3 (0) when /3 « 1. 



2.4. Computational units 

In our simulations, it is convenient to scale the fluid variables so that they are of order 
unity. To this end, we adopt the time unit 

t = IT 1 « 1.4M 8 rf yr, (20) 

where Mg = Mbh/(1O 8 M ) and r% = r/(10 3 r s ). As already noted, this is very short com- 
pared to the nuclear timescale of a massive star, t nuc = 0.007 M*c 2 / L Edd (M*) ~ 10 6 yr, and 
short even compared to the main-sequence Kelvi n-Helmholtz timescale, which approache s 



~ 3000 yr from above for very massive stars (e.g.. lBond et al.lll984l ; iGoodman &: Tanll2004l ). 
It is in part this disparity of timescales that causes us to doubt the ability of stellar feedback 
to stabilize the disk. More relevant to our simulations is the Kelvin-Helmholtz timescale for 
a radiation-pressure-dominated cloud of mass M* if we scale its radius R* by its Hill radius 



We choose our mass unit as 



3 V /2 ( k B \ _ in *J ( V® 



M °={^) (^|J= 10 ' 24M H,J • (22) 

Notice that, apart from the molecular weight /i, this is entirely determined by fundamental 
constants (/i Q ~ 0.62 is the molecular weight of a fully ionized gas at solar abundance). This 
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choice simplifies the Eddington relation between the mass of a homogeneous self-gravitating 
sphere and its gas-pressure fraction: 



2 



M^47M ^O(^1 . (23) 



P 2 \» 

At M = Mq, for example, this predicts (3 = 0.96. By no coincidence, M is characteristic of 
a moderately massive star. 

Because of the importance of self-gravity, it is convenient to choose the length unit so 
that Newton's constant is of order unity. The choice 

Lo ss ( 2 ^) V3 = 2.59 x 10" (^) V3 r, M r cm (24) 

implies G = (2tt)~ 1 . Then the units for surface density So, velocity vq, internal energy per 
unit area Uq and 2D pressure Pq are combinations of to, Mo, and L : 

S = M L 2 = 3.05 x 10 5 U^J r 3 - 2 M 8 - 4/3 g cm" 2 , (25) 
t; ^ = 5.75xl0 6 (^) r- 1/2 M~ 1/3 cm s" 1 . (26) 



\ t* J 

U = P = M to 2 = 1.0 x 10 19 ( — ^ r^ 3 M 8 " 2 erg cm" 2 . (27) 

Note that our length and time units, but not our mass unit, depend on radius. We 
might have avoided this by taking (K e M ) 1//2 and k 3 1 ^M 1 ^(27tG)~ 1 ^ 2 for our length and time 
units, respectively, whence our unit of surface density would be K~ l . At the radii of interest 
to us in a bright AGN disk, however, S ~ So ^> ft" 1 [eq. fl25|) ]. so Lq is the more convenient 
length standard. Also, Q enters more than once into the Euler equation ( |T2|) . A symptom 
of this is that the Hill radius works out rather simply: 

AL \ 1/3 fM^ 1/3 



R "^{u^) r ^ 376 {wJ u <- (28 > 

This is approximately the largest size at which a bound fragment of mass M* can withstand 
the tidal field. 

With these units, we express our dynamical variables in dimensionless form: 

£ = U=W> P = ( 2 9) 

S Uq Pq Vq 



- 11 - 



2.5. Accretion Rate 



As a consequence of the shearing-sheet boundary conditions, the joint average over time 
and space of the radial mass flux Hv x can be shown to vanish. Thus, we cannot expect 
to measure the mass accretion rate (M) directly in our simulations. Nevertheless, we can 
measure M indirectly from energy or momentum balance. 



At large ra dii in a stead y keplerian thin disk, energy balance is expressed by crT e 



cff 



3MQ 2 /8tt (e.g.. iPringld 1 1 9 8 lh . In our simulations, the effective temperature is a function of 
local parameters via the cooling function ( !T6|) . Therefore, one local estimator for M is 



M A 



47rA 



(30) 



On the other hand, steady angular momentum balance requires MQr 2 = T — Tq, where 
T is the "viscous" torque, and To is a constant that can be neglected at large radii. In a 
thin keplerian disk, this reduces to M = 37r^S, where v is the effective viscosity. In our 
self-gravitating shearing sheets, the role of T is played by r (G xy + H xy ), in which G xy and 
H xy are the offdiagona l compon e nts of the vertically integrated gravita tional and Reynolds 



37rz/S if we defind_| v = aP/VtY, 



stresses as defined by iGammid (120011 ) and iJohnson &: Gammid (120031 ). This agrees with 



2 G xy + H. 



a =3 



■'■fj 



and 



27i(G xy + H. 



xyj 



n 



(31) 
(32) 



In those simulations that reach a statistical steady state, the spat iotemporal ay erages 
of Ma and M a agree, and consequently {t c ) Xj y tt = (ar 1 )^^ -1 (e.g., IPringld Il98lh . We 
prefer the estimator ( )30|) in these steady cases because it fluctuates less than M a . When 
the sheet fragments into isolated masses that secularly cool, however, thermal equilibrium 
does not hold. Then M a is the more reliable estimator, and the "dissipation" of the mean 
shear associated with the stresses in eq. ( 1311) is balanced by increasing epicyclic motions of 
the fragments. 



1 As we use a different equation of state, our normalization of a differs from that of I Johnson fc Gammie 



(|2003|) by a factor involving the adiabatic index. 
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Numerical method 



Our simulations use a modified form of the code developed by iGammid (120011 ) and 
ohnson fc Gammid (120031) . This is a self-gravitating hydrodynamic code based on ZEUS 



( jStone fc Normanlll992l ). which is a time-explicit, operator-split, fi nite-difference m ethod on 



a staggered mesh. Details and tests of this code are described by IGammid (120011 ). We just 
emphasize some important points here. 



The code uses the standard "shearing box" boundary conditions (e.g., lHawley et al. 



19951 ). For a rectangular domain of dimensions L x x L y , all fluid variables / satisfy 
f(x,y,t) = f(x,y + L y ,t), 



f(x, y, t) = f(x + L x ,y- ptL x , t) , (33) 



except that v y (x,y,t) = v y (x + L x ,y 



-^VttL x , 



t) 



Poisson's equation is solved by 



discrete Fourier transforms. Mass is conserved up to round-off error, so the areal average of 



(Masset 


2000; 


Gammie 


2001) 



However, even without cooling, total energy — the sum of kinetic, internal, gravitational, 
and tidal energy (the tidal potential 4>t = \^lx 2 ) is not conserved. Part of the reason is the 
shearing-periodic boundary condition ( 133]) . which maintains the mean shear. The nearest 
thing to an energy integral is the Jacobi-like quantity 



d 3 xES(x) 



U_ 



(34) 



but this is not const ant unless cooli ng exactly balances the dissipation of mechanical energy 
by the stresses (13"T1) (IGammid 1200 ll ). 



Furthermore, there are numerical errors. ZEUS's energy equation is not written in flux- 
conservation form; indeed, such a form may not be possible for a razor-thin sheet whose self- 
gravity is described by the three-dimensional Poisson operator (1T41) . The velocities and mass 
densities are updated in such a way that the changes in kinetic energy due to gravitational 
forces are slightly inconsistent with changes in the self-gravity itself. The errors are second- 
order in space but only first-order in time. When density fluctuations are well resolved by 
the grid, these errors are minor: typically less than 10% over several hundred dynamical 
times. For compact fragments that span only a few cells, however, the fractional error in the 
binding energy after several shearing times can become more than 100%. Simulations of Jeans 
collapse in nonshearing (Q — > 0) sheets exhibit the problem clearly when the Jeans length is 
comparable to the grid scale. Massive (M ^> Mo) fragments are especially problematic with 
our "soft" equation of state, because thermal and gravitational energies nearly ca ncel for a 
nonrotating bound object that is radiation-pressure dominated (f3 <C 1). Because iGammie 



13 



( 120011 ) and iJohnson fc Gammid (120031 ) use a relatively "hard" equation of state, the error is 
less important, although it may contribute to the small deviations from the expected relation 
in their Figure 12. In our simulations, however, we found that the error could cause spurious 
dissolution of (originally) bound fragments. 

We therefore corrected the error as follows. At every time step, we calculate the expected 
change in the Jacobi integral (134)) due to cooling and to the stresses G xy + H xy acting on 
the mean shear. This is accurate to first order in the time step At. The change in T 
computed by ZEUS over the same step is slightly different. The error can be predicted to 
O(At) in terms of the state variables at the beginning of the step and the finite-difference 
algorithm that updates them. We compensate for this predicted error by multiplying the 
internal energy in every cell by a common factor^ This is equivalent to an extra cooling 
or heating term. The required change in the internal energy is never more than 1% per 
time step in the simulations reported here. Based on simulations of Jeans collapse and other 
tests, we believe that this procedure is sufficient to identify the boundary in parameter space 
between sheets that fragment and those that do not. Unfortunately however, the residual 
energy errors prevent us from following the merging of fragments to very large masses. 
Lagrangian methods, such as N-body methods and Smooth Particle Hydrodynamics, avoid 
this particular source of error because the non-dissipative parts of the numerical equations 
of motion, in cluding; gravitational term s, are fundamentally hamiltonian and have energy 



integrals (e.g. iMonaghan fc Pried l200ll ). 



When fragmentation occurs, the local cooling time t c = U(t) /A can become very short 
in the low-density regions between fragments, which contain very little mass. To prevent 
rapid cooling from limiting the time step, the internal energy is updated according to the 
stable scheme 

u(t + At) = YT~A«7i? < 35 > 

so that that U remains positive for all At. 



4. Results 

In this section, we present results from our simulations. The black-hole mass is taken 
to be 10 8 M Q , except in §4.41 where M B h = 10 9 M Q . §4.11 and §4.21 illlustrate nonfragmenting 



2 It might be better to use a local correction based on the density gradient and velocities, but since the 
gravitational interactions are intrinsically nonlocal, we were unable to find a satisfactory measure of the local 
error. 
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and fragmenting regimes, respectively. In §4.3[ we characterize the boundary between these 
regimes more systematically. Except where otherwise noted, all simulations are performed 
for domain sizes L x = L y = 10L at resolution N x = N y = 256. Experimentation indicates 
that varying these numerical parameters upward or downward by factors ~ 2 makes little 
difference to the results, at least through the early stages of fragmentation. 

We start all simulations with uniform surface density and internal energy, but with 
small random velocity perturbations added to the equilibrium velocity field v = —1.5Qxe y . 
Apart from these perturbations and from the numerical parameters cited above, the initial 
conditions are therefore characterized by three physical parameters in addition to the black 
hole mass: the reference radius ro, or equivalently the angular velocity f2 = (GMbk/vq) 1 ^ 2 ; 
the initial surface density, and the initial internal energy per unit area, t/j. 



Figure [T] shows the evolution of several diagnostic quantities in a simulation for ro = 
10 3 r s , Ej = 2S , and U = lOt/o- The quantities shown in the plots are implicitly averaged 
over the grid, and in some cases weighted by mass. When it is important to be explicit, we 
use angle brackets with appropriate subscripts, e.g. 



Occasionally overbars are used instead of brackets when we want to emphasize the time 
dependence of a spatial average, and we indicate by context or in some other way whether 
weighting by mass has been applied. Double brackets ((...)) a or ((. . .)) M will indicate 
averages over time as well as space. 

In this notation, the areal average = (X)a = £j since mass is conserved by our 
equations. For these choices of ro and £ mentioned above, we find that a statistical steady 
state is reached after t ~ 200f2~ 1 in which all of the quantities shown in Figure [T] fluctuate 
around long-term averages ((. . .}) that appear to be independent of the initial value U, as 
long as U is not so low that the sheet immediately fragments without cooling. In particular, 
((Q))a ~ ({Q))m ~ 1- Unless otherwise noted, we prefer to start from a hot state Qi ^> 1. 

To verify the steady state, we compare the cooling and heating rates. Averaged over 
the time interval between 200 and 800 the cooling time ((t c )) ~ 29.9f2 _1 , and 
the viscosity parameter ((a)) ~ 0.035, so that ((t c ))Q « ((a))^ 1 as required by thermal 
equilibrium ( §2.5)) . The accretion rate [Panel (b)] is somewhat larger than the Eddington 
rate Mgdd = O.22M 8 e _1 M yr _1 for this black hole for canonical values of the global radiative 



4.1. Case I: No permanent fragments 
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Fig. 1. — Time evolution of a shearing sheet at mean radius ro = 10 3 r s about a 10 8 M Q black 
hole, with mean surface density (S)a = 2X , box size 10L x 10L , and resolution 256 2 . 
(See §2.41 for definitions of the units L and S .) Panel (a), solid lines: Averaged Toomre 
parameters (Q)a (black); (Q)m (red). Dashed lines: averaged gas-pressure fractions {13) a 
(black); (f3) M (red). Panel (b): Average accretion rate (M)a via eq. (1301) [M yr -1 ]. Panel 
(c): Average cooling time (£ c )a (red) and maximum surface density T, max (black). Panel (d): 
Average internal energy per unit area (U)a (black) mid-plane temperature (T(0))a (red). 
No permanent fragments form, so areal and mass-weighted averages (. . .)a.m are similar. 
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efficency e = Ld^/Mc 2 « 10%. However, because r enters the shearing-sheet equations 
only via Q, these results could be mapped to r ~ 200rs around a 10 9 M Q black hole, where 
the accretion rate would be sub-Eddington. Panel (c) shows occasional strong peaks in 
surface density; however (Q) , which is a reciprocal measure of midplane density, never drops 
far below unity, showing that no permanent fragments form. Notice that the mass-weighted 
average {Q)m is systematically less than the areal average (Q)a because p(0) correlates with 
E. 



4.2. Case II: Permanent fragments 



In a constant-M, constant-a disk, Q declines with increasing radius (e.g. iGoodman 



20031 ). making fragmentation more likely. In a simulation with (E)a = 0.7So at tq = 
4 x 10 3 r s , the disk cools to Q < 1, and fragments form; as the sheet passes through {Q)a = 1 
the dimensionless cooling time Q(t c )A ~ 0.5. Starting from a uniform state, the mass 
first concentrates into azimuthal filaments, which then fragment into several dense clouds. 
After merging, a single bound object containing most of the mass results from our standard 
10Lo x 10Lo simulation (Fig. [2]). Unable to collide with itself, the object steadily cools. It 
shows no tendency to subfragment, as its Kelvin- Helmholtz time (12~T1) is longer than Q^ 1 , 
which in turn is longer than its internal dynamical time. 

No steady state results since we omit fusion reactions (which would ignite below the 
resolution of our grid). For the same (S)a/So, however, fragmention is avoided at the 
slightly larger radius 3 x 10 3 r S) where we measure ((M)) ps 2M yr _1 , congruent with an 
Eddington-limited disk feeding a 1O 8 M black hole at 10% radiative efficiency 

Figure [2] displays the fragmenting simulation at t — 53f2 _1 , after the dominant fragment 
has coalesced. As shown by the lower left panel, Q < 1 within the fragment, meaning that 
its midplane density is well above the Roche value. About half of the rest of the sheet is 
also dense at the midplane, but these regions have very little mass. The mass in the bound 
object is about 60M w 614M , 86% of the total. At this time, {{Q)) M ~ 0.019, and 
((/?)) m ~ 0.45. Note that this j3 is larger than what we expect for a nonrotating Eddington 
model of the same mass (eq. 123]) . This may in part be a numerical effect of our finite spatial 
resolution, which softens the gravitational force: as the radial density profile in Fig. [3]) shows, 
the half-mass radius of the object is approximately two cell widths. Another cause of the 
discrepancy may be the large rotational kinetic energy of this object, T/|W| ~ 0.13. 

The energy of this object is partitioned as follows: thermal energy = 936£'o; kinetic 
energy (measured with respect to its center of mass, mainly rotational) E k = 672E ; tidal 
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Fig. 2. — A bound fragment at t — 53f2 _1 in a simulation for r = 4x 10 3 r s and (S)^ = 0.7S . 
Clockwise from upper left: Surface density £/E ; potential vorticity £ (£ = 0.71 in the initial 
uniform state); gas pressure fraction (3; and local Toomre parameter Q [eq. (BJ]. 
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Fig. 3. — Black curve: Density profile of the fragment in Figure [2] versus distance r c from its 
center. Some 99.9% of the bound mass lies at r c < 0.46Lo, whereas the Hill radius for 6OOM 
is ~ 1.46 L . Solid red curve: Radial profile an Eddington model of mass 60.1M and outer 
radius 0.46L projected into two dimensions. Dashed red curve: A projected Eddington 
model of the same mass having the same central surface density as the fragment. The object 
appears to be distended by a combination of rotation and numerically softened gravity. 
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potential energy E t = —0.8E ; gravitational self-energy E p = —50Q0E . Thus the total 
energy of the object in its center-of-mass frame is negative, implying that the object is 
bound and doomed to contract indefinitely. 

Our 2D approximation facilitates merging because fragments cannot avoid one another 
vertically. The importance of this can be judged by examining the epicyclic motions of frag- 
ments, if one is willing to assume that the vertical and horizontal epicyclic amplitudes scale 
together. The epicyclic energy per unit mass of a free particle in the Keplerian shearing sheet 
is x 2 + (2y + 3Qx) 2 is a constant. For an isolated fragment, the corresponding characteristic 
quantity is 

E cpi = X -M [vl + {2v y + 3{lx) 2 ] . (37) 

Here M is the mass of the fragment, while the overbars mark the position and velocity of 
its center-of-mass. We define the (radial) epicyclic amplitude by R cpi = ^2E cpi / (MVL 2 ). 
The importance of the third dimension for collisions can be judged by comparing R cpi to the 
physical radius of the object, R*, or to its Hill radius, Rh [eq. f l28|) ]. We presume that the 
latter is the more relevant comparison, at least until R*/Rh < 0.1, because objects within 
one another's Hill sphere undergo a complicated motion in 3D that allows many opportunities 
for close passage. 

At t = 48. 4f2~ 1 in the above- described simulation for r = 4 x 10 3 r s , there are two frag- 
ments, with masses 50.6Mo and 2Mq, so that the Hill radius associated with their combined 
masses is Rh ~ 1.4 Lq. The epicyclic amplitude for the smaller mass is i? ep i ~ 3.14Lo- Later 
in our 2D simulation, these two fragments merge. We conclude that the merger might have 
been delayed or perhaps even avoided in 3D. 

In order to explore merging among more fragments, we have performed a simulation 
for the same ro, (£)aj and resolution as in Figure [2] but with four times the standard box 
size, i.e. L x = L y = 40Lq and NX = NY = 1024. The first bound fragments appear at 
t ~ Along one filament, fourteen small fragments form, and eight merge in pairs. It is 

easily shown the two-body problem decomposes in the shearing sheet into uncoupled motions 
of the center-of-mass and relative coordinates, as in free space. Therefore, adopting the 
approximation that the epicyclic motions of the two components of each pair are uncorrelated 
until shortly before they merge, we add their epicyclic amplitudes in quadrature, r cpi = 
( r e P i,i + r epi,2) 1,/2 5 an d compare this to the Hill radius based on their combined mass, i?n = 
[G(mi + m2)/3f2 2 ] 1 / 3 , with the results shown in Tabled] The data in the last column suggest 
that these encounters might have proceeded somewhat differently in three dimensions. 
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(m 1 + m 2 )/M m 2 /mi r cpi /R u 



2.17 


0.61 


3.10 


3.17 


0.83 


2.40 


1.37 


0.99 


1.15 


1.43 


0.96 


1.68 



Table 1: Epicyclic amplitudes of merging pairs. 
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Fig. 4. — Derived properties of statistically steady, nonfragmenting simulations versus radius 
and scaled surface density [E; eqs. (|25|) & ( 129|) ]. for M B h = 1O 8 M . The various symbol 
types mark corresponding simulations in all four panels. The cases for £ = 1,2 fragment 
beyond 2 x 10 3 rs & 3 x 10 3 rg, respectively, hence are not shown. The averages of Q and /3 
[panels (a) &(b)] are mass- weighted and systematically smaller than the corresponding areal 
averages. 
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Fig. 5. — The fragmentation boundary. Panel (a): Red dots mark simulations that frag- 
mented, black dots mark those that did not. Panel (b): M (solid line) and a (dashed line) 
along the fragmentation boundary, i.e. for the uppermost black dots in Panel (a). These 
represent maximal rates of gravitational transport without fragmentation. For 0i & t Cj i, see 
the text. 




Fig. 6. — Comparison of M along the fragmentation boundary with G03. Dashed curve: As 
in Panel (b) of Fig. (j5J). Solid curves: Predicted values of M for the indicated Q & a taken 
from the simulations, using eqs. (11) & (A3) of G03. 
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4.3. General picture of the self-gravitating regime 

In this section, we summarize the general trends found in our shearing-sheet models, 
with particular attention to the conditions for fragmentation. 

As discussed above, the eventual statistical steady state of our simulations, when it 
exists, are defined by two control parameters, mean surface density, (X), and angular ve- 
locity, f2 = (GMbh/^o) 1 ^ 2 - We have explored the parameter ranges 0.5S < (£) < 10S 
and 10 3 rs < r$ < 5 x 10 3 , with Mbh = 10 8 M Q . Throughout most of this regime, the ef- 
fective Shakura-Sunyaev parameter is measured to be a > 1CT 2 , so that self-gravity would 
likely have dominated the angular momentum transport even had we included MHD in our 
simulations. The corresponding accretion rates are 0.01 < M < 20 M & yr -1 . 

FigureHJshows several steady-state quantities as functions of our two control parameters. 
With the physics in our models, steady states are not possible after fragmentation, so these 
quantities are measured from simulations that did not fragment. As seen in the first panel, 
the mass- weighted average of Q is typically slightly less than < 1. Recall that our definition 
of Q [eq. 01])] is simply a reciprocal measure of the midplane density relative to the Roche 
density; under adiabatic compression by a nonlinear density wave, the internal energy of 
the gas rises in step with its density, so that a bound fragment may be avoided even if Q 
falls briefly below unity. Mass weighting tends to emphasize these transiently compressed 
regions. The areal average of Q is typically > 1. At a given £ = (E/E ), Q decreases 
toward larger radii, where cooling is stronger. Panel (b) shows that the mass-weighted /3 
increases with increasing radius and decreasing surface density. Panel (c) shows that M 
is much more sensitive to surface density (£) than to radius or, equivalently, to Q. For 
comparison, the Eddington rate for a black hole of 10 8 M 8 M Q is MEdd ~ 2M 8 eo.iM yr -1 , 
where e = O.leo.i = L^^/Mc 2 is the global radiative efficiency of the disk. For £ = 1, the 
Eddington rate is achieved at tq rs 3 x 10 3 r s . At the same M, however, our models fragment 
when r >4x 10 3 r s . The final panel shows that generally a increases with increasing radius 
at fixed S. But a has a complicated dependence on surface density. The largest value 
encountered in any of our non-fragmenting simulations was a max ~ 0.4. 

Figure |5] shows the boundary between those cases that fragment and those that do not 
in the plane of (r , X), our two control parameters. Each dot represents a simulation, all done 
with L x = L y = 10Z/Q and NX = NY = 256. We have checked that the boundary between 
the fragmenting (red) and nonfragmenting (black) cases is not significantly altered at higher 
resolution (NX = NY = 512). Higher resolution does make a difference, however, when the 
most unstable wavelength is short, which happens when the accretion rate (and thus surface 
density) is very small: our standard resolution begins to fail at M < O.1M yr" 1 . 
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Panel (a) of Figure [5] shows that the maximum dimensionless surface density that the 
disk can support without fragmenting declines rapidly with increasing radius. Over the 
range 10 3 < r /rs < 4 x 10 3 , the boundary can be fit to a power law: Ef rag ~ 6(r/10 3 rg) _L5 . 
Since, as shown in Figure HI the local accretion rate is much more sensitive to E than to r , 
it follows that M max also declines swiftly with radius; this is confirmed by the second panel. 
In fact, M max declines some two orders of magnitude between 10 3 rs and 5 x 10 3 rs- 

The gas pressure fraction 0i and cooling time t Cj i marked in the second panel are neither 
mass-weighted nor areal averages: instead, they are computed for the same r and (E) 
as those in the simulation, but with Q set to unity in the equation of state rather than 
its measured steady-state value. This allows us to compare the observed fragmentation 
boundary with Gammie's criterion (Qt c ) crit = constant = 0(1). Except at the innermost 
radius shown, we find that the cooling time is indeed O^" 1 ) along the boundary. However, 
fragmentation occurs at r = 10 3 rs when (E) > 5 even though Qt c ^> 1; this is possible 
because <C 1, so that bound fragments are only marginally stable against collapse even 
without energy loss. For Mbh = 10 8 M Q , this regime is reached only at local accretion rates 
far above the Eddington rate, but not so for larger Mbh, as will be shown in §4.41 

We have compared our simulations with the alpha-disk models of G03 for the case that 
viscosity is proportional to total pressure. Figure 1 of G03 displays curves of constant a and 
Q in a plane of M versus r. For Q = 1 and plausible a, there are generally two branches 
to the curve: a high-M solution, which has high surface density and low 0, and a low-M 
solution, which has the oppositie properties. These branches join at r m 10 3 rg, so that 
Q > 1 for all solutions at smaller radii. To compare with these predictions, we take the 
measured values of a and mass-weighted Q from the simulations along the fragmentation 
boundary shown in Panel (b) of Fig. [5j and we insert these values into the model for M 
from G03. The results are shown by solid lines in Fig. There are again two solutions for 
M at each radius, with radiation pressure dominating the upper (higher M) solution, and 
gas pressure dominating the lower. But while Q is roughly constant along these curves, a 
is not — a decreases rapidly with decreasing radius. The actual M directly measured in the 
simulations (dashed curve) lies above the upper branch at r < 4 x 10 3 rg and has slightly 
higher surface density. The differences between the predicted and measured values of M 
may be due in part to the assumption of uniform conditions in the a models, so that mass 
and areal averages differ. 
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4-3.1. Scaling to other black-hole masses 

Apart from fundamental constants and numerical parameters (grid resolution, domain 
size, etc.), the statistical steady states of our self-gravitating shearing sheets are entirely 
determined by just two control parameters^] f2 and (£). Therefore, although we have fixed 
Mbh = 10 8 M Q in our simulations and studied the outcomes as functions of r and (£), we 
can scale our results to other black-hole masses by recasting them in terms of the two control 
parameters above. For ease of writing, we will omit the angle brackets from (£) and (£) 
henceforth. 



Our most important result is the fragmentation boundary. As noted above, £f rag pa 6r 3 
with x ~ 1.5. Since r§ oc Mbh and £q oc f2 4//3 [eq. fl25|) ]. this can be recast as 



, N -(2+x) 
r Z/6 



2xl0 6 [M^r 3 ) gem" 2 . (38) 



^frag 

We compare this with the surface density required for accretion at the Eddington rate. 
From Panel (c) of Fig. m it appears that M is much more sensitive to £ than to radius. This 
implies that £ oc r -2 at fixed M. We will attempt to explain this scaling below, but for the 
moment, we simply accept it. Since M increases by a factor 10 2 as £ increases by 4, we 
estimate M oc T, y with y ps 3.3. Figure H] also indicates that £ ~ 0.7 yields M p^ 2M Q yr~ 1 , 
which is the Eddington rate for Mbh — 10 8 M Q and radiative efficiency e = 0.1. This 
coincides with £f rag at r 3 ps 4 for M 8 = 1. Rewriting the relation M ps 2 M yr _1 (£/0.7) y 
as M oc (£f2 4//3 ) y , we find that the radius beyond which a self-gravitating accretion disk will 
fragment if it accretes at a fraction fn of the Eddington rate is, taking x — 1.5 and y = 3.3, 

r crit w 4 x 10 3 M 8 ~°- 87 m a2 r s w 0.04 M 8 °- 13 m a2 pc. (39) 



For comparision, iGoodmanl ( 120031 ) 's equation (10) predicts that Q = 1 in an alpha disk at 
r 3 p^ 2.7(a/M 8 ) 2 / 9 if the viscous stress is proportional to total pressure and <§C 1. This is 
roughly half of eq. ( 13"9~|) for M 8 = 1 and a = 0.4 (the largest value found in our simulations), 
but the scaling with black-hole mass is different. As Panel (b) of Figure S] shows, however, 
is closer to 1 than to at r cr i t for M 8 = 1, so precise agreement is not to be expected. 

We promised to discuss why £ oc r~ 2 at fixed M and Mbh- When self-gravity controls 
the accretion rate, Q ~ 1, so that the midplane density p(0) oc r~ 3 . It follows from vertical 



3 Actually, the metallicity of the gas should be counted as a third parameter. Since it enters the opacity 
(Tl8|) as well as our mass unit (|22]l . it cannot be entirely scaled out of the simulations, for which we have 
taken fi = Hq throughout. 
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radiative and hydrostatic equilibrium that the half thickness of the disk is 

^«^(l + 2 g- 1 )- 1 (l- j 9)- 1 (40) 
one 

to the extent that (3 is vertically constant. This gives the familiar result that h ~ constant 
in steady disks where radiation pressure dominates. Then we would have E = 2hp(0) oc r -3 , 
not r~ 2 , for constant Q. Fig. HJshows, however, that (3 > 0.7 at r 3 > 1 for E = 0.7, the value 
that gives a roughly Eddington accretion rate for Ms = 1 . Th us h may vary with radius 



through the factor (1 — 0) 1 . Now eq. (A3) of iGoodmanl (120031 ) predicts that 

-9/10 



(1 - /3 )-i / 3 1 /2+(6-i)/io w .35(a .iM 8 )- 1 /io rh -4/5 ( _^ ) ,.21/20 ? (41) 



where b = or b = 1 according as i/ oc P or v oc P ra d: clearly it makes little difference 
to the value of (1 — 0) when this is < 0.5. Although a is not constant with radius in our 
self-gravitating models, the dependence in eq. (I4TI) is so weak that (1 — [3)~ l and hence h 
are approximately linear in r when (3 > 0.5. This explains why E oc r -2 , but it also shows 
that this scaling holds only over a limited range of r and Mbh- 



Equation ( 1391) shows that self-gravity is important at a smaller multiple of r§ for larger 
Mbh at a given Eddington fraction m; it then follows from eq. (]4ip that the self-gravitating 
regime is characterized by smaller /3 for larger Mbh- In fact, for Mbh ^ 10 9 M , we estimate 
that < 0.1 at r cr i t , so that fragmentation may occur with little cooling, as demonstrated 
in §4.41 Thus, while it remains true that the local dynamics of a self-gravitating disk is 
determined by Q and E, the particular scaling (13U1) . which depends upon our power-law fit 
to the fragmentation boundary over a limited range of is likely to be modified for black- 
hole masses much above 10 8 M . On the other hand, for black holes much less massive than 
our fiducial value, Kramer's opacity will dominate over electron scattering at r cr i t ; in view 
of the sensitivity of the radiation fraction (14"T|) to k, this also will modify eq. ( 13"9"|) . Thus, the 
latter equation is probably quantitatively reliable within only a narrow range around M 8 = 1. 
Nevertheless, the trend is surely correct: namely, that r cr it, the radius beyond which accretion 
at the Eddington rate would cause fragmentation, occurs at a smaller multiple of rg for larger 
M BH - 



4.4. Fragmentation at long cooling times 

As noted in the softening influence of radiation pressure on the equation of state 
may allow fragmentation even for Qt c ^> 1. One way to enter this regime is to increase the 
mean surface density. At fixed radius and fixed Q, the gas pressure fraction (3 w 0.5E~ 3 / 2 



- 26 - 




(b) m = 53 



(c) m = 60 



Fig. 7.— Simulation for r = 621r s , (S) A = 9.6S of 1O 9 M SMBH. Line plots: Evolution 
history of accretion rate calculated from eq. (130]) . a, and mass-weighted Q and (3. Because 
the disk is not in thermal balance, the M shown here does not reflect that due to turbulent 
stress, which would be smaller. Lower left: Three fragments form at time t = 53f2 _1 . Lower 
right: A single fragment is left after mergers at time t = QOQ^ 1 . 
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[eq. IA8] , while the cooling time t c oc kS 2 . For example, with Q — 1 and E = 6 at ro = 10 3 r s 
and M BH = 1O 8 M , our equation of state yields (3 = 0.03 and fit c = 244. At r < 200r s 
along the fragmentation boundary shown in Fig. 0, (3 is already very small and t c Q is large. 
The surface density exceeds what is required to support the Eddington luminosity but might 
occur if mass were dumped into the disk by a violent event such as a merger. 

Another way to enter this regime is by increasing the mass of the black hole at a 
fixed Eddington fraction, m. At Mbh = 1O 9 M and rh = 1, the disk is still strongly 
dominated by radiation pressure at the smallest radii where self-gravity is important. Setting 
Z E /eo.i = rh = 1, M 8 = 10, and a = 0.01 in equations (A2)-(A4) of G03, we estimate that 
Q = 1 occurs at r = 621r s , where (3 = 0.02 and E = 9.6. 

We have done one simulation with these values of Mbh, r o, and E (Figure [7]). The 
simulation starts at Q = 1.2 and cooling time t c Q = 49. The large but declining values of M 
in the first panel are computed from the thermal equation (1301) rather than the stress equation 
( 132|) . which would predict M m in the initial phases when a ~ 0. A sustained balance 
between radiative cooling and turbulent heating is never achieved, even though the cooling 
is slow and proceeds smoothly to Q < 1. After roughly one cooling time, at t ~ 50f2 _1 , the 
sheet collapses to an azimuthal filament. This quickly breaks into three massive fragments, 
which merge into one at 60f2 _1 . The energies of the final object are = 4.6 x lO 5 ^, 
-Ekin = 7.5 x 10 4 i?o, E ti d = — 15.1-Eo, and E grav = —5.4 x 10 5 E Q , respectively, so that it is 
marginally bound. 

As this example shows, the mass scale for fragmentation is very large if it occurs at 
<C 1. This is to be expected from the Eddington quartic (1231 . The fragments inherit 
an initial (3 similar to that of the disk because their formation occurs roughly adiabatically 
at high specific entropy, T 3 /p ps constant. The final object weighs 893M ~ 9000M Q at 
(f3) m ~ 0.07, close to the prediction from eq. ( 12"3"|) but larger than the value /?disk ~ 0.02 
expected for a uniform sheet at Q — 1 with this mean surface density. 

We have not attempted a thorough a survey of parameter space for Mbh = 1O 9 M as we 
did for 1O 8 M . One obstacle is that the cooling time becomes very long, especially at small 
radii. Another is that the gravitational stress remains small up to the point of fragmentation, 
in contrast to the situation for 10 8 Mbh where a max ~ 0.4 (Fig. |5]). At a < 10~ 2 , heating 
by MRI becomes important, which cannot be explored with this 2D code. MRI might have 
prevented fragmentation for the disk parameters of Fig. [7J since Q > 1 is predicted for these 
parameters if a > 10~ 2 , and fragmentation did not begin until Q < 0.5. The point, however, 
is that self-gravity alone was not able to supply a sufficiently large a despite the long cooling 
time. Thus our simulations demonstrate that fragmentation can occur at dimensionless 
cooling times Qt c ^> 1 when radiation pressure dominates, ^ < 1. 
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Discussion and Conclusions 



Radiation pressure and appropriate opacities are part of the minimal physics needed to 
study gravitational turbulence and fragmentation in near-Eddington AGN accretion disks. 
We have included these and used them to test the predictions of Goodman (2003) for the 
maximum radius at which AGN disks can support steady accretion at the Eddington rate 
without fragmentation. We are in qualitative agreement with that paper, but quantitatively 
we find that the critical radius is about twice as large as was predicted for Mbh = 1O 8 M . 
We are also generally in agreement with lGammid ( 120011 ) 's criterion for fragmentation, except 
that fragmentation may occur for Qt coo \ ^> 1 if radiation pressure dominates. 

Beyond this, however, our local 2D approximation limits what we can explore and what 
we can conclude. Magnetohydrodynamic (MHD) and magnetorotational (MRI) processes 
cannot be represented, at least not directly. This probably does not much alter the conditions 
for marginal fragmentation, because the effective viscosity due to self-gravity is much larger 
in this regime than what MRI can provide. We cannot, however, exclude the possibility 
that MRI, or some other effective viscosity that does not respect the conservation of specific 
vorticity, might interact with the self-gravity in subtle ways, for example by promoting 
secular instabilities at Q > 1, or by enabling transitions among states of different M at the 
same £ and Q; we saw evidence for the latter behavior when we added an ad-hoc viscosity 
v oc P gas to our code, but we have not been able to understand those results and therefore 
have not presented them here. Our approximations also cannot represent magnetized winds 
or global spiral arms, which might in principle remove angular momentum at rates enhanced 
by ~ r/h compared to transport within the disk, allo wing; a lower and hence less self- 



gravi t ating surface density for th e same accretion rate ( iGoodmanl 120031 ; iThompson et al. 
20051 : iHopkins fc QuataertlboiOah . 



Nor can we test lGoodman &: Tan! (120041 ) 's suggestion that fragments grow rapidly up to 
the isolation mass, which is typically ~ 10 5 M Q . There are two principal obstacles. One is 
numerical: very large self-gravitating masses are very strongly radiation-pressure dominated, 
and therefore only marginally bound when in virial equilibrium, so that small energy errors 
can cause spurious expansion or collapse. This is likely to be a difficulty for many numerical 
algorithms besides ZEUS in low-/?, self-gravitating regimes. The second is physical: our 
2D results show that in shearing sheets where multiple bound fragments are present, the 
gravitational interactions between fragments quickly increases their epicyclic motions to 
amplitudes exceeding their Hill radii , so that they would be expected to scatter into the 
third dimension if that were allowed ( iRafikov &: Slepianll2010l ). 



Notwithstanding these limitations of 2D, our results strongly suggest that if the disks 
of bright QSOs extend at constant M to > 0.01-0.1 pc, then bound objects will form with 
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individual masses of at least 3OOM , and possibly much more. Collectively, these "stars" 
will dominate the local surface density of the disk, though they may be accompanied by an 
optically thick layer of distributed gas. The stars will attain epicyclic dispersions bounded 
by Safronov numbers = GM^/(R^ { ) < 1, so that even if they contract to their main- 



sequence ra dii and are stable enoug h to remain there, (M*/Mbh)(?V-R*) <C 1 unless M* 



> 



10 5 M Q (see iGoodman fc Tanl 12004 for a review of the nominal main-sequence properties of 



very massive stars). Physical collisions will be important unless or until the objects collapse 
to blac k holes. One is thus lead to imagine a model for the disk similar to that advanced lon g 
ago by ISpitzer fc Saslawl (Il966[ ) , and more recently by iMiralda-Escude fc Kollmeierl ( 120061 ) . 
We imagine formation of (very massive) stars within a disk, however, rather than formation 
of a disk from a pre-existing dense nuclear star cluster. 

Are there any observational signatures of such a fragmented disk that might distinguish 
it from the conventionally imagined smooth one? One such signature may be the super-solar 
metallicity inferred from the broad lines, which appears not to correlate with the ge neral 
star formation rate in the host traced by far-infrared emission (ISimon fc Hamannll2010l ). and 
therefore may implicate formation within the nuclear disk itself. Another signature might be 
deviations from the spectral energy distribution expected from a steadily accreting, optically 
thick disk. IGoodman fc Tanl ( 120041 ) pointed out that the viscous accretion time at r ~ 10 3 rs 
is typically somewhat less than the minimum main-sequence lifetime, so that massive stars 
formed there might — if they are sufficiently stable — migrate to the inner edge of the disk 
before dying. In that case, if these stars dominate the surface density and are not fully 
enshrouded by diffuse gas, the local color temperature of the disk might be intermediate 
between that of the stars themselves and the effective temperature implied by the accretion 
rate. Gravitational microlensing observations are beginning to test the variation of apparent 
disk size with wavelength on relevant lengthscales; the evidence is consistent with color 
temperature oc r~ 3 / 4 as expected for steady, opti cally thick disks, but suggests that the d isks 
are larger at a given wavelength than expected (jPooley et al.l 120071 ; iMorgan et al.l 120101 ) . 
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A. Effective 2D Equation of State for vertically constant f3 

Vertical hydrostatic equilibrium in the combined gravitational fields of the central mass 
and of the disk itself is described by 

1^ = -n 2 z - AttG f p(z')dz' , (Al) 
pdz Jo 

Putting p = K((3)p^ 3 [eq. Q] and adopting the dimensionless Emden variables 

p(z)=p(0)9\ £ = z/h, h 2 = nG * Q)2/3 , (A2) 

leads to 

ie +e - _ 4^(o) = _ ?• (A3) 

Using the initial conditions 0(0) = 1 and 0'(0) = 0, eq. fl A3 j) can be reduced to a quadrature: 

i = 2^\^ dW (A4) 

io V2Q + 1 + (1 - w 2 ) + (1 - W 2 ) 2 + (1 - w; 2 ) 3 

The height-integrated density and pressure become 

S = 4v / 2p(0)/i/ 3 (Q), P = 4v / 2A'(/3)p(0) 4 / 3 /i/ 4 (Q), (A5) 

in which 

4(Q) - Z 1 (1 ~ W2) ^ = 3, 4. (A6) 

7o V2Q + 1 + (1 - w 2 ) + (1 - w 2 ) 2 + (1 - w 2 ) 3 

For all w G [0, 1] and Q > 0, the denominators of the elliptic integrals flA6j) vary by 
at most a factor of 2. Hence we approximate these integrals with single-point Gaussian 
quadrature scheme, 

(l-w 2 ) k f(w 2 )dk*u k (wl), (A7) 



o 



in which the point Wk e [0, 1] and the weight Uk > are chosen so as to make this scheme 
exact for functions f(w 2 ) = A + Bw 2 with arbitrary constants A and B. For fc = 3 and 
k = 4, the two integrals are close enough that the same value of Wk can be used for both; 
this leads to the approximations (jHJ), which are accurate to < 1% for all Q > 0. In practice, 
we prepare a table with a certain range of E and £/, within which we calculate the integrals 
iE} directly. For conditions outside the table, the code uses the approximations (jSJ). 
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Eliminating h between eqs. flA2j) and the first of eqs. (IA5I) and then expressing p(0) in 
terms of Q via eq. (jlj) leads to eq. (JBJ). Using this to eliminate K and /i from the second of 
eqs. flA5p then yields equation (J3J) for P in terms of Q and S. But P is related to the internal 
energy by eq. (jSJ), so (jSJ) can be recast as (j7|). Finally, after replacing K with its explicit form 
03]), equations ([6]) and (J7]) can be rewritten in terms of the dimensionless variables introduced 
in §21] as 

2 15 (l~/3) _^Q* 128 [I 3 (Q)] 3 fr 

/? 4 ~[/ 3 (Q)] 6 ' ^ 3Q/ 4 (g) S3" (Abj 

Equations (1A8I) implicitly determine /3 and Q given [/ and S, as exemplified by Fig. [HI after 
which P follows from eqs. (J^J) or (JSJ). 



REFERENCES 

Bartko, H., et al. 2010, ApJ, 708, 834 

Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825 
Collin, S., & Zahn, J. 1999, A&A, 344, 433 

Davies, R. I., Sanchez, F. M., Genzel, R., Tacconi, L. J., Hicks, E. K. S., Friedrich, S., & 
Sternberg, A. 2007, ApJ, 671, 1388 

Dhanda, N., Baldwin, J. A., Bentz, M. C, & Osmer, P. S. 2007, ApJ, 658, 804 

Gammie, C. F. 2001, ApJ, 553, 174 

Ghez, A. M., et al. 2003, ApJ, 586, L127 

Goodman, J. 2003, MNRAS, 339, 937, (G03) 

Goodman, J., & Tan, J. C. 2004, ApJ, 608, 108 

Hamann, F., & Ferland, G. 1993, ApJ, 418, 11 

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 

Hirose, S., Blaes, O., & Krolik, J. H. 2009, ApJ, 704, 781 

Hopkins, P. F., & Quataert, E. 2010a, ArXiv e-prints 

-. 2010b, MNRAS, 1085 



-32 - 



010" 




(a) 



10" 



10 



(b) 



_E=0.01L 

- z ~ z -o 

_E=100L„ 



10 10 

u/u 



10" 



10" 



10" 









_U=0.01U Q 












U=100U 



10"' 



10" E/2L o 10 



10' 



(c) 



=3-10 




10 



10"' 



10 u/u 10 



(d) 



10' 



Fig. 8. — Solutions to equation (1A8|) for (3 and Q in the parameter space S and U . In panel 
(a) and (c), we show the behavior of Q and (3 as a function of S for a fixed value of f/. 
Different lines are for different values of U as shown in the plots. In panel (b) and (d), we 
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different lines are also shown in the plots. The four plots are the equation of state that we 
use in our code. 
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